subroutine Possion2DT3CR1
    use Mesh2DT3CR_m
    use FES2DT3CR1_m
    use SMatrix_m
    use IterSSolver_m
    implicit none
    type(mesh2DT3) :: mesh2D0
    type(mesh2DT3CR) :: mesh2D
    type(FES2DT3CR1) :: fs
    type(SMatrix) :: A
    type(IterSSolver) :: solver
    real*8, allocatable :: u(:), RHS(:)
    character*80 :: fileName
    integer, parameter :: label(4) = (/1, 2, 3, 4/)
    integer :: iter
    
    call mesh2D0.square(5, 5, 3)
    call mesh2D.initMesh(mesh2D0)
    call fs.initFES(mesh2D)
    call A.initSMatrix(fs%nDof)
    allocate(u(fs%nDof), RHS(fs%nDof))
    RHS = 0.0D0; u = 0.0D0

    call fs.assembleHeatStiff(1.0D0, COMMON_2D, A)
!    call sortSMatrix(A)
!    call displaySMatrix(A)

    call fs.assembleHeatSource(1.0D0, COMMON_2D, RHS)
!   write(*, *)"RHS"
!   write(*, *) RHS
    call fs.applyBCZero(A, RHS, label)
    write(*, *) "starting computing..."
    call solver.initIterSSolver()
    call solver.PCGSSOR(A, RHS, u, 1.5D0, 1)

    write(*, *) "solution: "
    write(*, *) "min: ", minval(u), " max: ", maxval(u)

    fileName = "CRPossion.dat"
    call mesh2D.exportSolution(fileName, u)

    call A.deallocateSMatrix()
    deallocate(u, RHS)
end subroutine
